################################################################################################
#### 运行配置文件
################################################################################################

source /public/home/xxf2019/20231121_singleMuti/config/config.sh
conda activate archr_v2

################################################################################################
#### 结合RNA的细胞分类，计算atac按照501bp计算peak
################################################################################################
#### call peak
${Rscript_archr} ${scripts_dir}/run_archR.callPeak.R \
--rna_data_file ${input_dir}/testis_combined.Rdata \
--input_dir ${input_dir} \
--cluster_cell_file ${config_path}/cluster_celltype.csv \
--cpu 10 \
--out_path ${output_dir}/atac_res 

################################################################################################
#### 标记并展示RNA聚类结果
################################################################################################
## 注释细胞类型
${Rscript_archr} ${scripts_dir}/annotation_celltype.R \
--input_file ${input_dir}/testis_combined.Rdata \
--cell_type_file ${output_dir}/atac_res/testis_merge_all_qc_barcode_new.txt \
--out_path ${input_dir} 

## 展示细胞亚群及marker基因
${Rscript_archr} ${scripts_dir}/rna_celltype_show.R \
--input_file ${input_dir}/testis_combined.annotationCellType.Rdata \
--out_path ${output_dir}/rna_celltype 

################################################################################################
## 质控ATAC细胞并重新call peak
## RNA的异常细胞也同时去掉
################################################################################################
## 去除异常的一簇细胞，该Patchytene只在testis01患者中有比例为11%（700多个细胞），远高于别的患者的2.5%（300多个），去除以后比例正常
${Rscript_archr} ${scripts_dir}/atac_cell_check.R \
--comine_data_file ${output_dir}/atac_res/testis_combined_peak.combineRNA.Rdata \
--rna_file ${input_dir}/testis_combined.annotationCellType.Rdata \
--out_path ${output_dir}/qc_atac
## ${output_dir}/qc_atac/testis_combined_peak.combineRNA.qc.Rdata
## ${output_dir}/qc_atac/testis_combined.annotationCellType.qc.Rdata
## 20852个细胞

## 继续质控，生殖细胞和体细胞分布重聚类，去除离群的细胞
## 体细胞和生殖细胞分开
${Rscript_archr} ${scripts_dir}/atac_cell_check.v2.R \
--comine_data_file ${output_dir}/qc_atac/testis_combined_peak.combineRNA.qc.Rdata \
--rna_file ${input_dir}/qc_atac/testis_combined.annotationCellType.qc.Rdata \
--out_path ${output_dir}/qc_atac_v2
## 20606个细胞 = 10833胚系 + 9773体细胞
## ${output_dir}/qc_atac_v2/all/testis_combined_peak.combineRNA.qc.Rdata
## ${output_dir}/qc_atac_v2/all/testis_combined.annotationCellType.qc.Rdata
## ${output_dir}/qc_atac_v2/germ/testis_combined_peak.combineRNA.qc.Rdata
## ${output_dir}/qc_atac_v2/germ/testis_combined.annotationCellType.qc.Rdata
## ${output_dir}/qc_atac_v2/somatic/testis_combined_peak.combineRNA.qc.Rdata
## ${output_dir}/qc_atac_v2/somatic/testis_combined.annotationCellType.qc.Rdata


################################################################################################
#### 分所有的、胚系和体细胞
## 重新call peak和鉴定motif
for type in ${type_list[@]}
do
echo $type

## call peak
${Rscript_archr} ${scripts_dir}/archr_callPeak.R \
--comine_data_file ${output_dir}/qc_atac_v2/${type}/testis_combined_peak.combineRNA.qc.Rdata \
--out_path ${output_dir}/qc_atac_v2/${type}

## 导出atac的矩阵,每个基因算同一类细胞的均值
${Rscript_archr} ${scripts_dir}/get_genescorematrix.R \
--comine_data_file ${output_dir}/qc_atac_v2/${type}/testis_combined_peak.combineRNA.qc.Rdata \
--out_path ${output_dir}/qc_atac_v2/${type}

## 导出RNA的表达矩阵,每个基因算同一类细胞的均值
${Rscript_archr} ${scripts_dir}/get_expression.R \
--rna_file ${output_dir}/qc_atac_v2/${type}/testis_combined.annotationCellType.qc.Rdata \
--out_path ${output_dir}/qc_atac_v2/${type}
done

################################################################################################
## 既往研究已报道的差异基因在我们细胞类型中的GSVA
for input_file in `ls ${output_dir}/qc_atac_v2/all | grep Gene | grep tsv`
do
echo ${input_file}
o1=`echo ${input_file} | sed 's/.tsv//'`

for gene_list_file in `ls ${config_path}/reportDiffGene | grep list`
do
echo ${gene_list_file}
o2=`echo ${gene_list_file} | sed 's/.list//'`

out_name=${o1}_${o2}
${Rscript_expressionTime} ${scripts_dir}/gsva_list.R \
--input_file ${output_dir}/qc_atac_v2/all/${input_file} \
--gene_list_file ${config_path}/reportDiffGene/${gene_list_file} \
--out_name ${out_name} \
--out_path ${output_dir}/celltype_plot/gsva
done
done

################################################################################################
#### 细胞整体图谱
################################################################################################
##########################################
## 展示细胞UMAP及marker基因
## rna和atac均展示
${Rscript_archr} ${scripts_dir}/rna_celltype_show.R \
--rna_file ${output_dir}/qc_atac_v2/all/testis_combined.annotationCellType.qc.Rdata \
--atac_file ${output_dir}/qc_atac_v2/all/testis_combined_peak.combineRNA.qc.Rdata \
--out_path ${output_dir}/celltype_plot

##########################################
## 气泡图
${Rscript_archr} ${scripts_dir}/rna_celltype_dotplot.R \
--rna_file ${output_dir}/qc_atac_v2/all/testis_combined.annotationCellType.qc.Rdata \
--atac_file ${output_dir}/qc_atac_v2/all/testis_combined_peak.combineRNA.qc.Rdata \
--scriptPath ${scripts_dir}/scScalpChromatin \
--out_path ${output_dir}/celltype_plot

##########################################
## 差异表达基因
## 会输出满足pct的所有基因以及满足foldchange的差异表达基因
<< EOF
cpu=10
## 基因至少在多少比例的细胞表达
pct_list=(0.01 0.1 0.25)
## log2foldchange用于
logfc_list=(0.1 0.5 1 1.5 2 2.5 3 )
for pct in ${pct_list[@]}
do
echo $pct
for logfc in ${logfc_list[@]}
do
echo $logfc
${Rscript_archr} ${scripts_dir}/diff_expression.R \
--rna_file ${output_dir}/qc_atac/testis_combined.annotationCellType.qc.Rdata \
--scriptPath ${scripts_dir}/scScalpChromatin \
--pct ${pct} \
--logfc ${logfc} \
--cpu ${cpu} \
--out_path ${output_dir}/celltype_plot/diff_expression
done
done
EOF

type=germ
pct=0.25
logfc=1
cpu=10

## 差异表达(生殖细胞里面)
${Rscript_archr} ${scripts_dir}/diff_expression.R \
--rna_file ${output_dir}/qc_atac_v2/${type}/testis_combined.annotationCellType.qc.Rdata \
--scriptPath ${scripts_dir}/scScalpChromatin \
--pct ${pct} \
--logfc ${logfc} \
--cpu ${cpu} \
--out_path ${output_dir}/celltype_plot/diff_expression/${type}
## 画图
${Rscript_archr} ${scripts_dir}/diff_expression_plot.R \
--gobp_file ${output_dir}/celltype_plot/diff_expression/${type}/one_vs_other.${type}.pct_${pct}.logfc_${logfc}.GOBP.tsv \
--scriptPath ${scripts_dir}/scScalpChromatin \
--pct ${pct} \
--logfc ${logfc} \
--out_path ${output_dir}/celltype_plot/diff_expression/${type}/plot

## 差异atac
${Rscript_archr} ${scripts_dir}/diff_atac.R \
--atac_file ${output_dir}/qc_atac_v2/${type}/testis_combined_peak.combineRNA.qc.Rdata \
--scriptPath ${scripts_dir}/scScalpChromatin \
--cpu ${cpu} \
--out_path ${output_dir}/celltype_plot/diff_atac/${type}

##########################################
## 细胞质控图
${Rscript_archr} ${scripts_dir}/cell_qc_plot.R \
--rna_file ${output_dir}/qc_atac_v2/all/testis_combined.annotationCellType.qc.Rdata \
--atac_file ${output_dir}//qc_atac_v2/all/testis_combined_peak.combineRNA.qc.Rdata \
--scriptPath ${scripts_dir}/scScalpChromatin \
--out_path ${output_dir}/celltype_plot/cell_qc

<<EOF
## 卞师兄质控参数
GCN2 <- subset(x = GCN,subset = nFeature_RNA>200&nFeature_RNA<4000) #10852
GCN3 <- subset(x = GCN2,subset = nCount_RNA>500&nCount_RNA<20000) #7347
GCN4 <- subset(x = GCN3,subset = nCount_ATAC>1000&nCount_ATAC<20000) #5868
GCN5 <- subset(x = GCN4,subset = nucleosome_signal < 4) #5868
GCN6 <- subset(x = GCN5,subset = TSS.enrichment > 2) #5860
GCN1 <- subset(x = GCN6,percent.mt<50) #5592

#### ATAC质控
## 目前最靠谱的质控标准是TSS富集的得分(评估ATAC-seq数据的信噪比)和唯一比对数(比对数如果不够，那么该细胞也没有分析的价值)
## https://www.archrproject.com/bookdown/per-cell-quality-control.html
## 在 ArchR 中，通过过滤细胞被识别为 TSS 富集分数大于 4 且独特核片段超过 1000 个的细胞

## 我们的数据分布
## TSS > 4 :TSS 富集评分指标背后的想法是，与其他基因组区域相比，ATAC-seq 数据在基因 TSS 区域普遍富集，这是由于与启动子结合的大型蛋白质复合物。
## nFragments below 1,000 独特核片段的数量（即未映射到线粒体 DNA）。很简单 - 具有很少可用碎片的细胞将无法提供足够的数据来做出有用的解释，因此应该被排除。
## FRIP:The fraction of reads in called peak regions, 默认不进行质控,我们样本中最低为0.148
## 质控标准TSS > 4, nFragments > 1,000

#### RNA质控
## 皮肤癌的NG文章的标准
## First, cells were removed if they had fewer than 200 genes expressed, 
## fewer than 1,000 unique sequenced reads (unique molecular identifiers) or 
## greater than 20% of counts corresponding to mitochondrial genes. 

## 我们的数据分布
## percent.mito [0,0.25]
## nCount_RNA [591,49943]:nCount_RNA is the total number of molecules detected within a cell
## nFeature_RNA [501.00,9915.00]:nFeature_RNA is the number of genes detected in each cell
## 质控标准：线粒体<0.25, 50000 > nCount_RNA > 500, 10000 > nFeature_RNA > 500
EOF

##########################################
## 展示差异的peak及其所在的motif
## 提取所有的peak及其所属类型、提取差异的peak其所属类型、提取所有存在motif的peak、提取所有motif的位置
for type in ${type_list[@]}
do
echo $type
${Rscript_archr} ${scripts_dir}/archr_markerMotif.R \
--comine_data_file ${output_dir}/qc_atac_v2/${type}/testis_combined_peak.combineRNA.qc.Rdata \
--scriptPath ${scripts_dir}/scScalpChromatin \
--out_path ${output_dir}/celltype_plot/diff_peak/${type}
done

##########################################
## 所有样本中peak-gene分层哪集簇、每簇对应的peak富集在哪些TF、每簇基因富集在哪些GO通路
for type in ${type_list[@]}
do
echo $type
${Rscript_archr} ${scripts_dir}/archr_peak2gene.R \
--comine_data_file ${output_dir}/qc_atac_v2/${type}/testis_combined_peak.combineRNA.qc.Rdata \
--scriptPath ${scripts_dir}/scScalpChromatin \
--out_path ${output_dir}/celltype_plot/peak2gene/${type}
done


################################################################################################
#### 拟时序
################################################################################################
<<EOF
## 基于RNA
## 基于表达数据建立拟时序
${Rscript_archr} ${scripts_dir}/run_monocle.R \
--input_file ${output_dir}/qc_atac/testis_combined.annotationCellType.qc.Rdata \
--out_path ${output_dir}/monocole

## 整合基因的atac、rna以及拟时序的结果为Rdata
${Rscript_archr} ${scripts_dir}/run_combine_rna-atac.addPseudotime.R \
--rna_data_file ${output_dir}/qc_atac/testis_combined.annotationCellType.qc.Rdata \ \
--rna_data_monocle_file ${output_dir}/monocole/testis.monocle.Rdata \
--atac_data_file ${output_dir}/qc_atac/testis_combined_peak.combineRNA.qc.Rdata\
--time_diff_gene_file ${output_dir}/monocole/pseudotime_ordergene.tsv \
--gff3_file ~/ref/GTF/gencode.v32.annotation.gff3 \
--gene_ensg_file ~/ref/GTF/gencode.v32.gene_ensg.tsv \
--out_path ${output_dir}/monocole

## rna拟时序画图
## RNA本身的以及，与时序相关基因的表达和峰开放
${Rscript_archr} ${scripts_dir}/run_monocle_plot.R \
--rna_data_monocle_file ${output_dir}/monocole/testis.monocle.Rdata \
--comine_data_file ${output_dir}/monocole/testis_combined_peak.combineRNA.Rdata \
--time_diff_gene_file ${output_dir}/monocole/pseudotime_ordergene.tsv \
--out_path ${output_dir}/monocole
EOF

################################
## 基于ATAC的时序分析
## 最后使用
${Rscript_archr} ${scripts_dir}/run_trajectory.R \
--comine_data_file ${output_dir}/qc_atac_v2/germ/testis_combined_peak.combineRNA.qc.Rdata \
--out_path ${output_dir}/trajectory


################################################################################################
## 分每个细胞亚群,分别鉴定motif以及peak2gene以及TF-gene的对应关系
################################################################################################
##########################################
## 鉴定motif以及peak2gene
## peak2gene默认是250kb内
## 生殖细胞
export type=germ
cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' | grep -E ${germ_cell} | awk -F, '{print $1}' | xargs -P 5 -i sh -c '
sh ${scripts_dir}/chromvar_p2gene.sh {} ${config_path} ${type}
'
## 体细胞
export type=somatic
cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' | grep -E ${somatic_cell} | awk -F, '{print $1}' | xargs -P 5 -i sh -c '
sh ${scripts_dir}/chromvar_p2gene.sh {} ${config_path} ${type}
'

## 同一细胞类型的peak在不同样本间的差异，检查是否存在差异，用于展示质控
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' `
do
cluster=`echo ${line} | awk -F, '{print $1}'`

${Rscript_archr} ${scripts_dir}/peak_qc_plot.R \
--comine_data_file ${output_dir}/subcell/${cluster}/${cluster}.combineRNA.motif_peak2gene.Rdata \
--scriptPath ${scripts_dir}/scScalpChromatin \
--out_path ${output_dir}/celltype_plot/cell_qc/peakMA
done

## 每个细胞类型分别输出motif
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
${Rscript_archr} ${scripts_dir}/chromvar_getmotif.R \
--comine_data_file ${output_dir}/subcell/${cluster}/${cluster}.combineRNA.motif_peak2gene.Rdata \
--out_path ${output_dir}/motif_position
done

##########################################
## 用新的参数鉴定TF-gene的可能关系
## 在生殖细胞里面跑
export cor=0.5
export maxDelta=0.75
export out_path=${output_dir}/tf_regulators_${cor}_${maxDelta}

## 生殖细胞
export type=germ
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' | grep -E ${germ_cell} `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
sh ${scripts_dir}/run_TFregulators.use.sh ${cluster} ${config_path} ${cor} ${maxDelta} ${out_path} ${type}
done

## 体细胞
export type=somatic
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' | grep -E ${somatic_cell} `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
sh ${scripts_dir}/run_TFregulators.use.sh ${cluster} ${config_path} ${cor} ${maxDelta} ${out_path} ${type}
done

## 链接germ和somotic到all
mkdir -p ${out_path}/all
ln -snf ${out_path}/germ/* ${out_path}/all
ln -snf ${out_path}/somatic/* ${out_path}/all


## 所有细胞类型中阳性的TF合并到一张表格
## 头行
export out_path=${output_dir}/tf_regulators_${cor}_${maxDelta}/all

cat ${out_path}/cluster0/TFregulatorPlots_All/*tsv | head -1 | \
awk -F'\t' '{OFS="\t"}{print "TF",$0,"cluster","cell_type"}' \
> ${out_path}/Positve_TF-Gene.tsv
## 具体数据行
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
cell_type=`echo ${line} | awk -F, '{print $2}'`
echo ${cell_type}
for file in `ls ${out_path}/${cluster}/TFregulatorPlots_All/ | grep tsv`
do
tf=`echo ${file} | awk -F'_' '{print $1}'`
cat ${out_path}/${cluster}/TFregulatorPlots_All/${file} | grep -v putative_targets | \
awk -F'\t' '{OFS="\t"}{print TF,$0,cluster,cell_type}' TF=${tf} cluster=${cluster} cell_type=${cell_type} \
>> ${out_path}/Positve_TF-Gene.tsv
done
done

## 提取存在调控的TF-gene
cat ${out_path}/Positve_TF-Gene.tsv | grep -E "YES|putative_targets" \
> ${out_path}/Positve_TF-Gene.onlyTargetGene.tsv

## 提取存在调控的TF-gene的基因列表
mkdir -p ${out_path}/gene_list
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
cat ${out_path}/Positve_TF-Gene.onlyTargetGene.tsv | grep -w ${cluster}| \
awk -F'\t' '{OFS="\t"}{print $1,$2}' | tr '\t' '\n' | sort -u \
> ${out_path}/gene_list/${cluster}.gene_interest.list
done

## 计算存在调控的TF-gene基因之间表达的相关性，用于确定TF-gene的共表达
## 内存容易溢出
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
if [ ! -f ${out_path}/expression_correlation/${cluster}_expressionCorrelation.tsv ]
then
sh ${scripts_dir}/gene_correlation.sh ${cluster} ${config_path} \
${out_path}/gene_list/${cluster}.gene_interest.list \
${out_path}/expression_correlation
fi
done

## 注释存在调控的TF-gene表达相关性
export cpu=20
sh ${scripts_dir}/tf_gene.AnnoationExpressionCor.sh ${config_path} ${cpu} ${out_path}





##########################################
## 鉴定TF-gene的可能关系
## 默认代码
## TF选择在所有细胞类型间存在差异的
## 先跑,结果比较结果后面会进行调整
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
echo " sh ${scripts_dir}/run_TFregulators.sh ${cluster} ${config_path} " | \
qsub -N ${cluster}_"run_TFregulators" -l nodes=1:ppn=10,mem=100gb,walltime=240:00:00 -q batch -d ${qsub_log_path}
done

## 所有细胞类型中阳性的TF合并到一张表格
## 头行
cat ${output_dir}/tf_regulators/cluster0/TFregulatorPlots_All/*tsv | head -1 | \
awk -F'\t' '{OFS="\t"}{print "TF",$0,"cluster","cell_type"}' \
> ${output_dir}/tf_regulators/Positve_TF-Gene.tsv
## 其它具体相关系数行
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
cell_type=`echo ${line} | awk -F, '{print $2}'`
echo ${cell_type}
for file in `ls ${output_dir}/tf_regulators/${cluster}/TFregulatorPlots_All/ | grep tsv`
do
tf=`echo ${file} | awk -F'_' '{print $1}'`
cat ${output_dir}/tf_regulators/${cluster}/TFregulatorPlots_All/${file} | grep -v putative_targets | \
awk -F'\t' '{OFS="\t"}{print TF,$0,cluster,cell_type}' TF=${tf} cluster=${cluster} cell_type=${cell_type} \
>> ${output_dir}/tf_regulators/Positve_TF-Gene.tsv
done
done

## 提取存在调控的TF-gene
cat ${output_dir}/tf_regulators/Positve_TF-Gene.tsv | grep -E "YES|putative_targets" \
> ${output_dir}/tf_regulators/Positve_TF-Gene.onlyTargetGene.tsv

## 提取存在调控的TF-gene的基因列表
mkdir -p ${output_dir}/tf_regulators/gene_list
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
cat ${output_dir}/tf_regulators/Positve_TF-Gene.onlyTargetGene.tsv | grep -w ${cluster}| \
awk -F'\t' '{OFS="\t"}{print $1,$2}' | tr '\t' '\n' | sort -u \
> ${output_dir}/tf_regulators/gene_list/${cluster}.gene_interest.list
done

## 计算存在调控的TF-gene之间表达的相关性
## 内存容易溢出
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
if [ ! -f ${output_dir}/tf_regulators/expression_correlation/${cluster}_expressionCorrelation.tsv ]
then
echo " sh ${scripts_dir}/gene_correlation.sh ${cluster} ${config_path} \
${output_dir}/tf_regulators/gene_list/${cluster}.gene_interest.list \
${output_dir}/tf_regulators/expression_correlation" | \
qsub -N ${cluster}_"gene_correlation" -l nodes=1:ppn=40,mem=250gb,walltime=240:00:00 -q batch -d ${qsub_log_path}
fi
done

## 注释存在调控的TF-gene表达相关性
export cpu=30
export out_path=${output_dir}/tf_regulators
echo " sh ${scripts_dir}/tf_gene.AnnoationExpressionCor.sh ${config_path} ${cpu} ${out_path} "| \
qsub -N "tf_gene.AnnoationExpressionCor" -l nodes=1:ppn=40,mem=250gb,walltime=240:00:00 -q batch -d ${qsub_log_path}

## 提取存在peak-gene关系的motif
${Rscript_archr} ${scripts_dir}/run_TFregulator_getMotif.R \
--postive_tf_file ${output_dir}/tf_regulators/Positve_TF-Gene.onlyTargetGene.annoExpCor.tsv \
--scriptPath ${scripts_dir}/scScalpChromatin \
--input_path ${output_dir}/subcell \
--out_path ${output_dir}/tf_regulators


##########################################
## 用新的参数鉴定TF-gene的可能关系
## 更宽松，能鉴定更多已报道的可靠的TF
cor=0.3
maxDelta=0.7
export out_path=${output_dir}/tf_regulators_${cor}_${maxDelta}

for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
echo " sh ${scripts_dir}/run_TFregulators.use.sh ${cluster} ${config_path} ${cor} ${maxDelta}" | \
qsub -N ${cluster}_"run_TFregulators" -l nodes=1:ppn=10,mem=100gb,walltime=240:00:00 -q batch -d ${qsub_log_path}
done

## 所有细胞类型中阳性的TF合并到一张表格
## 头行
cat ${out_path}/cluster0/TFregulatorPlots_All/*tsv | head -1 | \
awk -F'\t' '{OFS="\t"}{print "TF",$0,"cluster","cell_type"}' \
> ${out_path}/Positve_TF-Gene.tsv
## 具体数据行
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
cell_type=`echo ${line} | awk -F, '{print $2}'`
echo ${cell_type}
for file in `ls ${out_path}/${cluster}/TFregulatorPlots_All/ | grep tsv`
do
tf=`echo ${file} | awk -F'_' '{print $1}'`
cat ${out_path}/${cluster}/TFregulatorPlots_All/${file} | grep -v putative_targets | \
awk -F'\t' '{OFS="\t"}{print TF,$0,cluster,cell_type}' TF=${tf} cluster=${cluster} cell_type=${cell_type} \
>> ${out_path}/Positve_TF-Gene.tsv
done
done

## 提取存在调控的TF-gene
cat ${out_path}/Positve_TF-Gene.tsv | grep -E "YES|putative_targets" \
> ${out_path}/Positve_TF-Gene.onlyTargetGene.tsv

## 提取存在调控的TF-gene的基因列表
mkdir -p ${out_path}/gene_list
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
cat ${out_path}/Positve_TF-Gene.onlyTargetGene.tsv | grep -w ${cluster}| \
awk -F'\t' '{OFS="\t"}{print $1,$2}' | tr '\t' '\n' | sort -u \
> ${out_path}/gene_list/${cluster}.gene_interest.list
done

## 计算存在调控的TF-gene基因之间表达的相关性，用于确定TF-gene的共表达
## 内存容易溢出
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
if [ ! -f ${out_path}/expression_correlation/${cluster}_expressionCorrelation.tsv ]
then
echo " sh ${scripts_dir}/gene_correlation.sh ${cluster} ${config_path} \
${out_path}/gene_list/${cluster}.gene_interest.list \
${out_path}/expression_correlation" | \
qsub -N ${cluster}_"gene_correlation" -l nodes=1:ppn=40,mem=250gb,walltime=240:00:00 -q batch -d ${qsub_log_path}
fi
done

## 注释存在调控的TF-gene表达相关性
export cpu=20
echo " sh ${scripts_dir}/tf_gene.AnnoationExpressionCor.sh ${config_path} ${cpu} ${out_path} "| \
qsub -N "tf_gene.AnnoationExpressionCor" -l nodes=1:ppn=40,mem=250gb,walltime=240:00:00 -q smp -d ${qsub_log_path}



##########################################
## 所有已报道的关键的TF在motif水平、其本身开放程度以及表达水平进行染色
${Rscript_archr} ${scripts_dir}/plotEmbedding_TFregulators.R \
--rna_file ${output_dir}/qc_atac/testis_combined.annotationCellType.qc.Rdata \
--comine_data_all_file ${output_dir}/qc_atac/testis_combined_peak.combineRNA.qc.Rdata \
--report_tf_file ${config_path}/TF_upregulated_during_spermatogenesis.csv \
--out_path ${output_dir}/report_tf


################################################################################################
#### 标记并展示已知的重要基因在rna及atac的热图
################################################################################################
## 基因分为总的，TF以及nonTF
## 展示按照基因的表达聚类以及按照基因的固定顺序
${Rscript_archr} ${scripts_dir}/pheatmap_knowngene.R \
--rna_file ${input_dir}/testis_combined.annotationCellType.Rdata \
--comine_data_all_file ${output_dir}/qc_atac/testis_combined_peak.combineRNA.qc.Rdata \
--known_gene_file ${config_path}/SCOS_known_genes.csv \
--out_path ${output_dir}/pheatmap_knowngene


################################################################################################
## 共享数据
mkdir -p ${output_dir}/shareData/peak-gene
mkdir -p ${output_dir}/shareData/motif/cisbp
mkdir -p ${output_dir}/shareData/motif/jaspar

cd ${output_dir}/shareData

## 细胞类型注释
cp -rf ${config_path}/cluster_celltype.csv ${output_dir}/shareData/

## 转录组数据
ln -snf ${input_dir}/testis_combined.annotationCellType.Rdata ${output_dir}/shareData/

## atac包含peak的信息
ln -snf ${output_dir}/atac_res/testis_combined_peak.Rdata ${output_dir}/shareData/

## peak-gene
for file in `find ${output_dir}/cluster_all_result | grep correlation_magic.Rdata`
do
ln -snf ${file} ${output_dir}/shareData/peak-gene
done

## motif
for file in `find ${output_dir}/motif | grep rda | grep cisbp`
do
ln -snf ${file} ${output_dir}/shareData/motif/cisbp
done
for file in `find ${output_dir}/motif | grep rda | grep jaspar`
do
ln -snf ${file} ${output_dir}/shareData/motif/jaspar
done